source entropy.tcl

# collecting data 
time {
Charmm create entropyCalculator 
entropyCalculator setup {
   dir entropy
   cmd charmm.inp
   output charmm.out
   compressOutput no
}
entropyCalculator start
entropyCalculator do {
    load rtf ~/simulation/common/toph19_eef1.1.gouy.inp
    load para ~/simulation/common/param19_eef1.1.gouy.inp
    load psf 2mag.psf
    load coor 2mag.crd
}
Traj create A entropyCalculator { pept_n1.dcd pept_n2.dcd pept_n3.dcd pept_n4.dcd pept_n5.dcd}
entropyCalculator genIC
array set ic [entropyCalculator getIC ]
set nDihe [llength $ic(phi)]
#array set Residues [groupby ic {residue backbone}]

# list to contain dihedral data(in degrees). 
# Will be converted to nframe*nDihe two-dimensional array in C.
set DiheData {}
set i 0
A eachframe {
    refillIC 
    array set ic [ getIC ]
    for {set j 0} { $j < $nDihe } {incr j } {
       lappend DiheData [lindex $ic(phi) $j]
    }
    incr i
}
set nframe $i

}

# Divide dihedrals into groups by the residue it belongs to.
# Neglects the correlation between groups for calculation speed.
array set groups [groupby ic { residue }]

# calculate entropy 
time {
foreach group [array names groups] {
  #uncomment the next line to use histogram method to do entropy calculation.
  #set result($group) [MIE_hist $group 2 DiheData $nframe $nDihe 0.0 10.0 ]
  # the arguments frome left to right are:
  #'group': a list of dihedral angle index to be used in entropy calculation.
  #'level(2)': consider correlation up to order 2.Higher order calculation needs much more data point to converge.
  #'DiheData': the list that contains the raw data. 
  #'nframe': number of frames inside DiheData
  #'nDihe': number of dihedral angles inside DiheData(must be correct! entropy calculation function use this to find data).
  #'start(0.0)': the dihedral angle values starts from ?
  #'window width(10.0)': the winow width you use to construct histogram.
  set nlevel 3
  if {[llength $groups($group)] <= $nlevel} {set nlevel [llength $groups($group)]}
#  puts -nonewline "$group ( $groups($group)):"
  #uncomment the next line to use Nearest Neighbor method to do entropy calculation.
  set result($group) [MIE_NN  $groups($group) $nlevel DiheData $nframe $nDihe 10]
#  puts "$result($group)"
  # the arguments frome left to right are:
  #'group': a list of dihedral angle index to be used in entropy calculation.
  #'level(6)': consider correlation up to order 6.NN method converge faster than histogram but is time consuming.
  #'DiheData': the list that contains the raw data. 
  #'nframe': number of frames inside DiheData
  #'nDihe': number of dihedral angles inside DiheData(must be correct! entropy calculation function use this to find data).
  #'kth nearest neighbor': similar to window width in histogram method. When k>5, increasing k will not change the result too much(<0.00001%) but will increase the  consumed time.
}
}
# output entropy 
set s 0.0
foreach group [lsort -integer [array names groups]] {
  puts "$group ( $groups($group)): $result($group)"
  set s  [expr $s+$result($group)]
}
puts "Total S: $s"

